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Abstract —Accurate, fast and robust power flow (PF) algorithms are important in distribution systems (DSs) due to their 
great importance on management of DS applications. From short term to long term horizons, PF equations should be 
considered into the models. In this paper, a linear programming model is proposed for PF calculations of DSs. To develop 
the model, some incidence matrices are proposed for modeling DSs based on graph theory. Moreover, a new linear ZIP 
load model is proposed based on small variations of voltage phase angle in DSs. The linear ZIP load model along with 
incidence matrices, lead to developing a linear programming model for PF problem. The proposed PF algorithm may easily 
handle distributed generation, meshed topologies, and unbalanced networks. Since the developed model is linear, it could 
save in computational complexity and execution time of PF calculations. Furthermore, proposed linear PF constraints 
could be used in many system studies a structural PF model for DSs. Simulation cases are provided to show the accuracy, 
robustness, and performance of the proposed PF algorithm. 

Index Terms —distribution power flow, linear programming, ZIP load model, graph theory, distributed generation. 

I. INTRODUCTION 

The modern smart distribution systems (DSs) require analytical tools providing them decision support in real-time management 
of system performance. Among them, power flow (PF) calculation is the most important analysis in DSs. Many of computations 
like network reconfiguration, state estimation, and volt-VAR optimization are based on accurate PF results [1,2]. On the other 
hand, DSs are widely affected by the new paradigms caused by smart grid (SG) concepts [2]. DSs are going to utilize more meshed 
topology and to operate as an active network utilizing widespread distributed generations (DGs), with bidirectional flows of power. 
Moreover, PF models should be capable to efficiently handle topology changes imposed by network reconfigurations. Hence, it 
seems that robust and efficient PF methods are vital for reliable real-time management of modern DS. 

Nevertheless, almost all PF algorithms are nonlinear and require iterations based solutions. The main nonlinearity source is for 
voltage and current relationship of general ZIP load models. However, PF solutions should be obtained by an iteration procedure. 
This process is time-consuming and sometime may be diverged in case of ill-conditioning systems. Hence, this may restrict the 
reliable application of PF in DS management. In real-time analysis and optimization of DS, computation time and robustness 
become important factors and analytic solutions would be more convenient. An important motivation of this paper is to propose a 
fast and high-accurate PF algorithm originally developed for DSs based on linear programming approach. Hence, time-consuming 
procedures, such as Lower Upper (LU) factorization and forward/backward substitution, and ill-conditioned matrices, will not be 
concerned. All of these advantages lead to the reliable application of PF solution in real-time DS management. 

PF methods of DSs may be categorized into Newton-Raphson (NR) [1,2], backward/forward sweep (BFS) [3,4], and direct 
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methods [5,6]. Some characteristics of DSs like high RfX ratios of lines, mixed overhead-underground cable sections, phase 
unbalances, and weakly meshed networks, impose some difficulties in applying of conventional PF methods directly for DSs 
[1,2,5,6]. Usually, NR based PF methods suffer from the ill-conditioning problem and low convergence for DSs. Although some 
methods were proposed to avoid divergence and ill-conditioning problems [1], it is still a time-consuming process. PF methods 
based on impedance/admittance matrices, may not efficiently utilize special configuration of DS. 

The BFS [5,9], is one of the widely used methods in DSs, since it takes full advantage of weakly meshed topology of DS. In [3] 
a PF method based on BFS is presented in which, loads are modeled as constant impedance. Generally, PV nodes are considered 
by a post analysis. However, the convergence of this method deteriorates as the number of PV nodes increases and even maybe 
diverged [4,7]. Implicit Z-matrix method is another PF method in which loads are modeled as equivalent injections [10]. A PF 
algorithm is proposed in [11] where network loops are handled by breaking the ties creating some dummy buses. This method is 
further improved by using graph theory in [12], But it still needs a post-analysis for PV nodes. A direct method proposed is based 
on defining two so-called BIBC and BCBV matrices [5] [6]. Actually, these matrices calculate both backward and forward sweeps 
in a single iteration. However, the direct method is still an iteration based model and could not handle three phase transformers. 

There are some linear PF methods in the literature [2,7,8]. A linear PF approximated model for a balanced DS is proposed in [2] 
which provides an upper and lower bound for voltages of system buses. A linear approximation of ZIP loads is proposed in [7] 
which makes the PF equations linear. In [8], the Z-matrix method is approximated by linear equations with coefficients based on 
minimization problem. However, modeling PV nodes is still an important issue in these linear PF methods. 

The main idea behind this paper is to propose some PF constraints originally developed for DSs to be applied in optimization 
models. Common PF constraints are based on nonlinear NR model which have some disadvantages when applied in DSs. On the 
other hand, BFS and directs PF methods have iterative procedures which limit their application in the optimization model. Hence, 
in this paper, a novel PF model is proposed for DSs which provides some linear PF constraints. The proposed linear PF constraints 
could be easily employed in optimization models of DSs from operation to planning horizons. The PF model is then completed by 
a new linear ZIP load model to establish a linear programming model for the PF problem in DSs. In order to develop the model, 
some incidence matrices are employed for modeling of DS enable to model meshed and unbalanced networks. Based on the 
incidence matrices, PF equations are modeled as linear equations. The linear ZIP load model is established based on small variation 
of the voltage phase angle in DSs which makes load injections as a linear function of bus voltages. The proposed linear PF model 
leads to elimination of complex matrix procedures and time consuming iteration based solutions. Briefly, paper contributions may 
be expressed as follows: 

• The proposed PF method, takes full advantages of the special structure of DSs based on the graph theory. For radial and 
meshed DSs, one incidence matrix and two incidence matrices are required for network modeling, respectively. Hence, the 
minimum number of incidence matrices are employed compared to the other PF methods. 

• The topology of DS is only modeled by the incidence matrices and line characteristics and load parameters are modeled in 
separate matrices. Hence, any topology variation or change in characteristics of the network could be applied shortly into the 
model. This eliminates any need to recalculation of admittance or impedance matrices. Hence, it could widely reduce execution 
times for online application of the PF algorithm. This feature could be useful for reconfigurations studies in which DS operator 
should decide about the optimal status of tie-line switches in the network. 

• A new linear ZIP load model is proposed which makes PF algorithm linear. It utilizes small variations of voltage phase angle 
in DSs to linearize ZIP load model. 

• The proposed PF method can handle PV nodes directly into the model and no post sensitivity analysis or in-program iterative 
reactive power calculation is needed. Also, it may efficiently handle unbalanced networks. This is particularly useful for 
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studies like state estimation in which a full consideration of unbalances in the DS is needed. 

In the rest of the paper, in section II the basis for modeling of DSs is presented. In section III, the proposed ZIP load model is 
introduces. In section IV, modeling of DG units are presented. Proposed PF linear optimization model is presented in section V. 
Finally, some conclusions are provided in section VI. 

II. Modeling of Distribution System 

In this section, modeling of DSs is represented which is an extension of our previous work in [13,14] with the same principles. 
A. Unbalanced line model 

The loading of a distribution feeder is inherently unbalanced due to a large number of unbalanced loads, and un-transposed and 
multi-phases line segments. Hence, line segments should be modeled using the impedance matrix, Z abc (1). The Z abc includes series 
and mutual impedances of lines. Consequently, the voltage at sending end and charging current of lines connected to the bus might 
be calculated as in (2) and (3), respectively [13,14]. 
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z z 

ba—n bb—n bc—n 


Z Z , Z 

ca—n cb—n cc—n 


y abc _y abc _j_ ^ abc jj abc 


rabc,sh _ 


= (i/2)i; abc v; abc 


(1) 

( 2 ) 
(3) 


Note that U abc just shows the series current of line. Also, I abc - sh shows capacitance current of lines modeled at head buses in n 
model. The (l)-(3) provide a complete model for voltage and currents of an unbalanced system. For any phase failed to represent, 
the corresponding row and column in related matrices should be considered as null-entries. 



Fig. 1. A typical distribution network 
B. Incidence Matrices 

Consider a typical distribution network shown in Fig. 1 with an arbitrary tree graph shown with solid lines. Network lines for 
tree graph are called as "branches" with / index. The number of branches ( L ) always equals to N- 1 in a tree graph. Other lines 
which do not contribute in the tree graph are called as "tie-lines" with m index. The total number of lines including branches and 
tie-lines are equal to L+M. 

The loads and all shunt elements in each bus may be converted into an equivalent current injection [5], For tree graph, by direct 
application of the KCL, the relationship between branch currents and bus injection currents may be obtained using C matrix in (4a) 
[14]. The matrix C is called 'branch current to bus current incidence matrix'. The dimension of C is Lx(N- 1) and Cu equals to one 
if the current of line l is composed by injection of bus i, as shown in (4b). 


M=lc][i] 


(4a) 
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The relationship between voltage drops of buses and branch currents may also be obtained using transpose of the C [14]. The 
C r is of (N-l)xL dimension which il ,h element shows the existence of l ,h line in the route between i th bus and the swing. Hence, the 
C matrix could be solely used for complete molding of radial DSs. Let the voltage drop of bus i to the swing bus is represented as 
A Vi = Vi - Vi. Then, the matrix of all voltage drops of system buses could be represented in a matrix form as (5). 

[av] = [cf[z][(y] (5) 


Matrix Z is the diagonal matrix of branches impedances and is of LxL dimension. For unbalanced cases, diagonal elements of 
Z, should be replaced by the 3x3 matrix of three-phase line impedance shown in (1). 

Meshed topologies 

Some DS may be energized from both sides due to high voltage drop in long radial lines. In addition, distribution feeders may 
contain loops created by closing normally open tie-switches. Tie-lines affect the flow of voltage and current in the network. For 
each tie-line from bus i to bus j carrying current u m , injection currents of head buses will be changed. Changes in injections of head 
buses in (4a), could be handled using D matrix [14]. The D is 'branch current to tie-line current incidence matrix' and is of LxM 
dimension. New relation of branch currents is shown in (6a) using D matrix (6b). l/ m is the vector of tie-line currents. The number 
of bus injection remains constant in (6a), but M unknown variables of tie-line currents have been added compared to (4a). 


[t/HC][/] + [D][t/„] 




01-10 0 
01 - 10-1 


(6a) 

(6b) 


Each tie-line constitutes a fundamental loop with some tree branches which make a system of independent equations. Voltage 
drops in network loops could also be calculated using transpose of D as shown in (7) [14]. 

[0] = [Df[Z][C/] + [Z„][i7„] (7) 


Z) T is of MxL dimension. Also, Zm is a diagonal matrix of tie-lines impedances represented with m index in (7). M unknown 
variables in (6a) could be handled using M new equations of (7) for voltage drop in M tie-lines. Hence, network tie-lines could be 
handled only by another incidence matrix of D. It should be emphasized that bus voltage drops in (5) remain unchanged since 
routes from buses to the swing don not changed. Therefore, (5)-(7) could completely model PF equations for each DS. 

III. Linear ZIP Load Model 

A well-known load model is the ZIP model in which, each load has been considered consisting of three major parts: a constant 
power, a constant current and a constant impedance [15,16]. Mathematically, the load current may be determined based on the 
voltage of load bus and rating condition of the load as (8). 

i x ° = aY\ + + y(s*/v t )* (8) 


Calculation of /-part and p -part injections of the load currents are nonlinear equations. In this section, a linear ZIP load model is 
proposed based on small variations of the voltage phase angle. 
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A. z-part injection 

Injection of z-part is calculated using constant admittance. The load admittance and its phase angle are calculated in (9a) and 
(9b), respectively. Hence, real and imaginary parts of z-part injection are calculated in (9c) and (9d), respectively. 


Y =G +jB =S*=P 0 -jQ 0 

(9a) 

XY = 9 = tg' 1 {B/G) = tg _1 (-<2 0 /P 0 ) 

(9b) 

/ z ' re =P 0 v rc +e 0 v im 

(9c) 

i z ' im =-Q 0 v re +P 0 v im 

(9d) 

B. /-part injection 


In order to calculate i- part injection of load, fixed power factor criterion of the load should be satisfied. The fixed magnitude is 

determined based on rating condition. But, the phase angle of the load should be modified according to the phase angle of voltage 

and admittance (8). Almost all linear ZIP load models have assumed that both magnitude and phase angle of /-part current to be 

fixed [5,7,16]. However, it may not be a correct assumption based on fixed power factor criterion of the load. The magnitude of i- 

part injection is determined as (10a). Therefore, real and imaginary parts could be calculated as (10b) and (10c). 

I i= \So\ = ylPo+Qo 

(10a) 

i Ure =I'cos(S +9) 

(10b) 

i i,im = / 1 sin(c> + 9) 

(10c) 

In DSs, considering the voltage angle of the substation as a reference, the imaginary part of the voltage is smaller than the real 

part, by several times [15,16]. Hence, the phase angle and magnitude of the imaginary part of voltages of system buses lead to 

zero. This leads following triangular approximations in (11) to be applicable in DSs. 


c°s(4) = l; sin(£ ; ) = S j 

(11a) 

s,=tg- i (vr/vr)-vr/v?^r 

(lib) 

Utilizing cosines and sinus of the sum of two angles to calculate (10) and applying approximations of (11), imaginary and real 

parts may be calculated as (12a)- (12b), respectively. 


i '- rc = I ' (cos(d?) —v im sin(0)) 

(12a) 

i Um = / ‘ (v im cos {9) + sin(#)) 

(12b) 

C. p- part injection 


The p- part injection is a nonlinear function of voltage (8). Taylor series is a well-known method to linearize the nonlinear 

equations around any point of interest. The function of 1/v is analytical around unity and may be approximated with 2-v. Hence, 

the injection of p-part could be calculated as (13a). Using some simple calculations, real and imaginary parts of the p-part injection. 

could be represented as ( 13b)-( 13c), respectively. 


* P =(5 0 /v i )*=( P o -7Go)(2-v f )* 

(13a) 

i v ’ K =P 0 (2-v re )+<2 0 v im 

(13b) 

;r- im =P oV im -Q 0 ( 2 - V re ) 

(13c) 
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Collecting all expressions for real and imaginary parts of currents in (9), (12) and (13), the real and imaginary parts of current 
can be calculated as (15a)-( 15b) by k coefficients (15). 

i, l0 - re =kl°v; e +k 1 °V™ +k'° (14a) 

•Kim _ £lo y re + + ( 14b ) 

Coefficients of proposed ZIP load model are shown with ki,..,kt > in (15). The k coefficients are constant numbers which can be 
calculated in offline analysis. The model of (14) is linear which makes the PF equations linear. 

k l =(a-y)P 0 (15a) 

k 1 ={a-y)Q {) -pi' sin(0) (15b) 

& 3 =/31'cos(&) + 2yP 0 (15c) 

k 4 =(y-a)Q 0 (15d) 

k 5 = {a-y)P 0 + /?/‘sin(#) (15e) 

k 6 =pI'sin(0)-2yQ Q (15f) 

IV. Modeling of Distributed Generation 

Distributed generations are becoming more prevalent in the power industry. They include a broad variety of technologies such 
as wind turbines, photovoltaic systems, fuel cells, micro-turbines and different types of energy storage systems. DGs are more 
connected to the network in distribution levels. Hence, it is necessary to provide analytical tools for DGs in DSs. 

In general, the energy source of DG units may be categorized into determined and undetermined energy sources. Besides, they 
may be connected to the grid via synchronous generators, induction generators, and power electronic interfaces. Various energy 
sources combined with different energy converters result in various models. DG units could be modeled either as a PQ bus or PV 
bus in PF studies [6,17,18]. Some DG models are shown in Table I [6,19]. 


Table I- DG models for power flow studies 


Energy source 

Energy converter 

Model 

micro-turbine 

synchronous 

gen. 

PQ 

gas turbine 

synchronous 

gen. 

PV 

wind 

Induction gen. 

PV 

Fuel-cell 

Inverter 

PQ 

photovoltaic 

Inverter 

PQ 


For a DG modeled as PQ bus, injection of DG may be calculated using (14). The voltage of DG modeled as PV bus could be 
calculated as (16) in which approximations of (11) is applied. 

vT =|^ v |cos(4) + y|^''|sin(^.)^^'' + jVTS t (16) 

Injection of the DG unit could be determined as (17a). By some simple calculations, real and imaginary parts injections of DG 
could be represented as (17b)-( 17c), respectively. 
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The developed model in (16)-(17) may completely handle a DG unit modeled as a PV bus. Note that the second term in the 
numerator of (17b) makes the equation nonlinear. But, due to low amounts of phase angle in DSs, one may ignore the terms 
associated with 8 in (17b). 

V. Linear Programming Model for Power Flow 

A PF model based on linear programming is proposed here using the developed models. Although using an objective function 
(OF) in PF problem seems to be needless; but, notable application of optimization models in system studies makes it justifiable. 
The OF should be defined based on the type of the problem, while the proposed linear constraints could be used as efficient PF 
equations developed for DSs. PF constraints based on admittance matrix or Newton-Raphson models suffer from poor numerical 
instability and slow convergence rate for DSs. However, the proposed PF constraints are defined based on DSs and could be 
efficient in the convergence of the model. 

A. Injection current of a bus 

The bus injection current consists of three major parts; injections due to connected load modeled in (14), injections due to shunt 
capacitance of lines connected to bus based on k model presented in (3), and injections due to capacitor bank modeled as constant 
susceptance, B cp [15]. Collecting all injection terms, total injection current could be calculated as (18). 


i i =if + j (Bf+Bf)v i 

(18) 

The (18) is just valid for PQ load buses with specified active and reactive powers. For PV buses, the injection current may be 

determined using (17b) and (17c). 

B. Optimization model 

The complete proposed PF optimization model is represented in (19). It consists of a problem-dependent OF, exact constraints 

of PF equations, approximated constraints of load model, and of PV buses. The OF in (19a) is defined as the minimization of 

summation of bus voltage drops. However, it could be defined differently, since it does not affect the validity of the PF model. 

OF:min|x(v i rc -vr) + (vr-vr)J 

(19a) 

s.t. linear exact PF constraints: 

»r=ZQ/, re +ZA,x: 5 v/,v/ 

i m 

(19b) 

^ ,m =Zor+ZAx;\v/,v/ 

i m 

(19c) 

vr-v? =Z c .v {*,«r -xy?}, vi,vi 

i 

(19d) 

vf-vf =z ci [R^r + vivv/ 

i 

(19e) 

0 = { R K ~ *X m }+ Ry: - xy:, Vm 

l 

(19f) 

o=z D i\ R y + x,u? }+ r jc + x,x :, 

(19g) 


s.t. linear approximated load (PQ buses) constraints: 

if =k lV f+(k 2 —B- h -S, cp )v‘ m +k 3 , Vi gQ n (19h) 


if' 16 =(PfS-qf)/ 
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(19i) 


s.t. linear approximated PV buses constraints: 


, re = p. pv /V ( . pv , Vi eQ pv 

(19j) 

(P^'S-qf'' )yVj. pv , V/ e Q pv 

(19k) 

v* =V, pv , Vi £Q pv 

(191) 

v ™ =v pv S, Vi eQ pv 

(19m) 


The constraints of (19b)-(19g) represent PF equations of the network based on models developed in (5), (6a), and (7). Since no 
approximations were used in moldering DS, these constraints could exactly model PF equations. Hence, they could be used as an 
efficient alternative to conventional Newton-Raphson or admittance matrix models for PF equations. The constraints of (19b)- 
(19c) relate the line currents to bus injection and tie-lines currents based on (6a). The equations of (19d)-(19e) represent the voltage 
drop of system buses with respect to the voltage of the swing bus developed in (5). The constraints of (19f)-(19g) express voltage 
drop in loops created by tie-lines based on (7). For radial DSs, the (19f)-(19g) and the terms associated with D matrix in (19b)- 
(19c), should not be considered into the model. 

Injection currents of PQ buses are modeled in (19h)-(19i) based on (14) and (18). Since the exact injection of a ZIP load in (8) 
was linearized in (14), some approximations are involved in these constraints. However, one could use the exact ZIP load model 
which result in a nonlinear PF model. It could be solved iteratively or using a nonlinear optimization model. In simulations, errors 
associated with linearized models will be analyzed. Note that, injections of DGs modeled as PQ should also be modeled using 
(19h)-(19i). 

The constraints of (19j)-( 19m) are proposed for PV buses. The (19j)-(19k), represent the injection of PV buses based on (17). 
Besides, the (191)-(19m) define real and imaginary parts of voltage based on (16), respectively. Limits of reactive power in PV 
buses are ensured by a post-analysis. After convergence of the problem, the limitations of reactive power will be analyzed for PV 
buses. In case of any violation, the PV bus is turned to a PQ bus and the reactive injection power is adjusted to marginal values. 
Then, the PF model will be resolved. 

VI. Simulation Results 

The proposed models are validated by various simulations. Linear optimization models are simulated in GAMS [20] and typical 
nonlinear algorithms are modeled in MATLAB. The ZIP coefficients are assumed as a — 0.025, (> = 0.51, y = 0.465 [16]. All 
models are implemented in a personal PC equipped with 2.93 GHz processor and 4 GB RAM. 

A. Proposed Load Model 

In order to illustrate the performance of the proposed linear ZIP load model, it was implemented in a typical PF problem. Three 
different linear ZIP load models were compared with the exact nonlinear model as follows: 

• Model A: Nonlinear exact ZIP load model of (8). 

• Model B: The linear model of [7] in which just the p-part term is linearized and the phase angle of /-part is not modified. 

• Model C: The linear model of [16] in which the ZIP model is approximated by a ZI model with new coefficients. 

• Model D: The proposed load model which both nonlinear terms of ZIP model were approximated. 

The IEEE-18 bus test system was used for the comparison of load models [13]. The system is a radial distribution system, 
containing 16 busses at 12.5 kV, and 2 busses at 138 kV. The error in the calculation of voltage magnitude in different load models 
is shown in Fig. 2. The least error is for the proposed ZIP model. The error is typically less than 0.002 pu and 0.02 degree for 
voltage magnitude and phase angle, respectively. The ZI load model shows a relatively high error. This drawback may arise from 
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coefficients of ZI model calculated for fixed conditions. This may restrict the application of the model for changing conditions. 
Even though both models of B and D are linear, but the proposed ZIP model shows better accuracy. 



Bus 

Fig. 2. Calculation error of voltage in different ZIP load models 

The accuracy of the load model depends on the type of loads. In order to show the capability of the proposed ZIP model, the 
coefficients of [i and y were changed from 0 to 1 in steps of 0.2. In each combination of loads, the maximum magnitude error of 
voltage for load models of B and D was calculated. Calculation errors are shown in Fig. 3. As it could be seen, the calculation error 
of the proposed load model is less in all combination of loads. Also, since the constant current part of the load does not modify in 
load model of B, the error tends to increase with the increase in />. Hence, the proposed load model may be utilized as a reliable 
linear load model for DSs in which constant current loads are dominant. 


-3 



Fig. 3. Calculation error of voltage in different ZIP load models 

Accuracy of linear ZIP load models is proportional to voltage drop. The more voltage drop, the less accuracy of load model. 
The accuracy of the proposed ZIP model is also examined in various voltage drops. In Fig. 4, the maximum calculation error is 
shown. Since the term of T/v' was approximated around the unity, calculation error is increased with an increase in voltage drop. 
The model is accurate enough to 0.2 pu voltage drop. Hence, the calculation error of the model D may be less enough in typical 
operation of DS, since almost all distribution operators avoid a voltage drop of 10%. The proposed linear ZIP model may bring a 
good compromise between accuracy and performance. 



Fig. 4. Maximum calculation error versus maximum voltage drop 
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Proposed Power Flow Optimization Model 

The proposed PF optimization model was implemented in some test systems. PF outputs are compared with the direct PF method 
of [5]. Utilized PF methods are as: 

• Ml: proposed linear optimization model. 

• M2: direct PF method of [5]. 

1) Eight-bus unbalanced network 

An 8-bus unbalanced distribution feeder [21] shown in Fig. 5 is used to show the capability of the proposed PF model for 
unbalanced cases. Table II represents PF outputs for different methods. The highest error is in the order of 10' 4 shows the accuracy 
of the proposed method compared to the direct method. In another simulation, the proposed ZIP linear load model was utilized 
into the direct PF method of M2. Its results are as the same with results of Ml. Low enors are due to the approximations of linear 
ZIP load model and PF constraints are accurate as BIBC and BCBV matrices. In the last two rows of Table III, execution time and 
the number of iterations are represented. Execution times for the M2 and Ml are 4.37 ms and 0.671 ms, respectively. The proposed 
method is more than six times faster than the BIBC method for this case. Also, the solution is obtained with no iteration, whereas 
the M2 needs 6 iterations for a tolerance of 10" 6 . The proposed linear optimization model eliminates the need for the iterative 
solution. This may lead to a great saving in execution time to be implemented in real cases in near real-time applications. 
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Fig. 5. Eight-bus unbalanced test network 

Table II- Results of the 8-bus unbalanced system 


PF method 

M2 


Ml 

Load model 

ZIP 

linear ZIP 

Bus Phase 

Ivl 

5 

Ivl 

5 

a 

0.9977 

-0.0088 

0.9977 

-0.0088 

2 b 

0.9939 

-2.0903 

0.9939 

-2.0903 

c 

0.9734 

2.0968 

0.9735 

2.0968 

o b 

0.9820 

-2.0823 

0.9820 

-2.0824 

c 

0.9502 

2.0930 

0.9503 

2.0930 

4 c 

0.9481 

2.0931 

0.9482 

2.0931 

5 c 

0.9699 

2.0969 

0.9700 

2.0969 

6 c 

0.9676 

2.0970 

0.9677 

2.0970 

7 a 

0.9953 

-0.0096 

0.9953 

-0.0096 

8 b 

0.9787 

-0.0088 

0.9786 

-0.0088 

Running time (ms) 

4.37 


0.671 

Iterations 

6 


1 


2) IEEE 33-bus test feeder 

The PF methods of Ml and M2 are also implemented on the IEEE 33-bus test system [22], The test system consists of 32 
branches and 5 tie-switches. At first, tie-switches were ignored. The voltage profile using Ml is shown in Fig. 6. Results of applying 
M2 are almost the same as Ml and the highest error is 0.0002 pu for bus #18. The simulations are performed in 2.981 ms and 0.781 
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ms respectively for Ml and M2. It shows the performance of the proposed PF model for the non-iterative solution of the problem. 
The OF in the first case is equal to 1.5909 which shows an average 4.7% voltage drop in each bus. The least voltage is 0.918 pu 
for bus #18. Tie-switches was considered in another case whose results are also shown in Fig. 6. Improvement of voltage profile 
is obtained considering tie-switches. The OF is equal to 0.9286 which shows a decrease in the voltage drop across the network. 
The average voltage drop is reduced to 2.7%. 

The capability of the proposed PF method for handling DGs was investigated in another simulation. Two DG units worked as PV 
buses are placed in bus #18 and #33. They wish to adjust the voltage at 1 pu. The rated active powers are equal to 0.035 pu and 0.1 
pu, respectively. At first, tie-switches were ignored in simulations. Voltage profile is presented in Fig. 6. It could be seen that the 
proposed method can successfully adjust the voltage for PV nodes at 1 pu. Reactive power injections are equal to 0.049 pu and 
0.067 pu for DGs in bus #18 and #33, respectively. The OF is equal to 0.5109 which shows a decrease compared to the previous 
cases including no DG. It shows an average 1.5% voltage drop in each bus. Finally, tie-switches and DGs were both considered. 
It’s expected that voltage profile to be even better than previous cases. The voltage profile is shown in Fig. 6. It shows more 
decrease in voltage drop compared to the case with no tie-switches. The OF, in this case, is equal to 0.4215 which shows an average 
voltage drop of 1.2%. 



Fig. 6. Voltage profile of IEEE 33-bus system 

An important issue that may lead the PF algorithms to be diverged, is high R/X ratios. In order to show the robustness of the 
method Ml , simulations were repeated for different R/X ratios. In all simulation cases, the PF algorithm successfully converged. 
For high R/X ratios, the error increases since load model approximations become less accurate for higher voltage drops. Another 
problem which may diverge PF algorithm is ill-conditioned equations. It usually occurs when the network contains some very 
short lines or long lines. Some impedances are multiplied by ten and some others are divided by ten in another simulation. In all 
cases, the convergence of the proposed PF method is observed. 

3) IEEE 123-bus test feeder 

In order to show the capability of the model for handling large unbalanced networks, it was implemented on the modified IEEE 
123-bus network [22], Outputs of PF are shown in Fig. 7. The bus voltage, especially for phase ’a’, is somewhere below the desired 
level of system voltage. The OF, in this case, is equal to 11.3151. In the second case, a DG unit was placed at bus #95. The rated 
power of DG is equal to 0.25 pu working in PV mode. Voltage profile in the new case is also shown in Fig. 7. Improving the 
voltage profile as well as adjusting the voltage magnitude of bus #95 is obviously observed in Fig. 7. The OF is equal to 5.6948 
which shows a decrease in voltage drop in the network. Also, the DGs inject reactive power of 4.582 pu, 3.056 pu, and 4.115 pu 
into the network for each 'abc' phases, respectively. Running time is 8.32 ms and 1.15 ms for methods of Ml and M2, respectively. 
Fow running time with no-iteration procedure again reveal the performance of the proposed PF model. 


































X IT 1 Fig. 7. Voltage profile of IEEE 123-bus system with and without DGs 

11 i Other DS studies could utilize linear PF models as well. Volt-VAR optimization, expansion studies, harmonic PF, and filter 
i ^ 0 placement are of such complexed-model studies. For instance, reconfiguration is a mixed-integer nonlinear important problem 
I ^ 1 which should be done in near real-time. Heuristic methods suffer from an un-guaranteed global optimal solution. The proposed 
11V linear PF constraints could be used for utilizing speed and robustness of linear optimization models. A major advantage of linear 
I “1 A optimization models is reaching global optima in a polynomial acceptable time. 

VII. Conclusion 

In this paper, a linear programming approach was proposed for PF in DSs based on some proposed incidence matrices. Incidence 
matrices were defined based on graph theory. They could completely model DSs and could easily handle any meshed topology. 
Furthermore, based on small variation of voltage phase angle in DS, a linear ZIP load model was proposed which makes PF 
algorithm linear. Finally, a linear optimization model was proposed for PF calculations. 

The proposed PF model was implemented in some distribution test systems. The results showed the robustness, accuracy, and 
performance of the model. Time-consuming processes such FBS is not needed. Furthermore, topology-based incidence matrices 
lead to the elimination of time-consuming procedures, such as LU factorization and also ill-conditioned matrices. All of these 
advantages lead to the reliable application of proposed PF algorithm in DS analysis. It can provide decision support for real-time 
optimization in fundamental and harmonic analysis of DS problems such as state estimation, DS reconfiguration, Volt-VAR 
optimization, harmonic analysis, power quality studies, and reactive power procurement. 

Nomenclature 

In the paper, single parameters, single variables, and matrices of variables and parameters are shown with capital letters, small 
letters, and capital letters in brackets, respectively as follows: 


Sets 

U 

Network buses. 

l 

Network branches (tree graph) 

m 

Network tie-lines 

superscripts 

im 

The imaginary part of a complex number 

lo 

Load 

re 

The real part of a complex number 

sh 

Shunt capacitance of line in pi model 

z, i, p 

ZIP parts of voltage dependent load 

Parameters 












































Branch current to bus injection incidence 
matrix. 

Bus voltage to branch current incidence 
D 

matrix. 

Z = R+jX Line impedance (branches and tie-lines) 

P Active power. 

Q Reactive power. 

Y = G+jB Parallel admittance of system buses. 

k Load coefficients in proposed ZIP load model. 

a, f, y Coefficients of ZIP loads 
d The phase angle of the voltage 

0 The phase angle of load admittance 

Variables 

U Branch currents. 

V Bus voltages. 

I Bus injection currents. 
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